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ABSTRACT - We revisit a little known theorem due to Beltrami, through which the integra- 
o ' tion of the geodesic equations of a curved manifold is accomplished by a method which, even if 
<t^ inspired by the Hamilton- Jacobi method, is purely geometric. The application of this theorem to 
^ ' the Schwarzschild and Kerr metrics leads straightforwardly to the general solution of their geodesic 
equations. This way of dealing with the problem is, in our opinion, very much in keeping with the 
^ ! geometric spirit of general relativity. In fact, thanks to this theorem we can integrate the geodesic 
i—i ] equations by a geometrical method and then verify that the classical conservation laws follow from 
^ ; these equations. 
(N ; 
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1 Introduction 

In the preface to his classic textbook "The Analytical Foundations of Celestial Mechanics" [1], A. 
Wintner wrote in 1941 ". . . even the classical literature of the great century of Celestial Mechanics 
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appears to be saturated with rediscoveries. . . " He then illustrated this remark with some examples of 
results that had been rediscovered more than once, retracing their history back to the first discoverer. 
This is also true of the field of differential geometry and especially of its applications to general 
relativity (GR), where Malcolm MacCallum has used the particularly appropriate term "literature 
horizon" to describe the loss of familiarity with results that accompanies each new generation of 
relativists. We propose that an old theorem on geodesies due to Eugenio Beltrami (1835-1900) [2] , 
which as far as we know has not been quoted in texts on differential geometry since the publication of 
the classic books by L. Bianchi [3] and L.P. Eisenhart [4], deserves to be rediscovered. In his theorem 
Beltrami showed that the geodesic equations can be integrated through a method substantially 
analogous to the method of separation of variables used in the integration of the Hamilton- Jacobi 
equation. 

In GR we have some important spacetimes which are exact solutions of the Einstein equations 
and whose metric tensor components are known explicitly in a given system of coordinates. Starting 
from these components, one can write the geodesic equations 1 

d 2 x l . dx l dx k , , 

and then try to integrate them to determine the paths of test particles. The Schwarzschild spacetime 
(whose timelike geodesies can be used to calculate the advance of the perihelion of Mercury) and the 
Kerr metric (representing the gravitational field outside a rotating body or of a mathematical black 
hole) are two important examples whose geodesies can yield important physical results. 

Two methods are typically used to integrate the geodesic equations. Either one starts with the 
Lagrangian equations of motion (obtained from a Lagrangian C given by C = (dx l /ds) (dx k /ds), 
where for timelike geodesies s may be identified with the proper time) or with the corresponding 
Hamilton- Jacobi equations, in both cases representing a mechanical system governed only by a kinetic 
energy term, in which the effects of the gravitational field are represented by the curvature of the 
spacetime associated with the metric which determines this kinetic energy function. Now one is 
dealing with a mechanical system again instead of pure geometry. In our eyes, this approach seems 
to be a step backwards with respect to the spirit of GR. The motion of a test particle in a gravitational 
field is interpreted as the motion of a free particle in a curved spacetime which turns out to follow a 
geodesic. On the other hand, a completely "geometric" integration of the geodesic equations can be 
performed without referring to the equivalent point particle mechanical system. Once the geometric 
problem has been solved, the constants of integration can be interpreted as physical constants that 
are the first integrals of the motion in the classic approach. For its historical interest we reproduce in 
Sect. 2 a concise translation of the relevant section of Beltrami's article, only updating its notation. 



2 Beltrami's integration methods for geodesic equations 

We recall Beltrami's method for obtaining the solutions of geodesic equations [2, Vol. I, p. 366], an 
extension of the Hamilton- Jacobi integration method used for integrating the equations of motion in 
dynamics. This method originating with Beltrami in the later 1800s does not appear in recent books 
[7]; the last exposition of this method is reported in Eisenhart's classic book [4, p. 59]. 
Let us consider an n-dimensional semi-Riemannian manifold V n whose metric is represented by 

d s 2 = g ih dx % dx h (i, h — 1, 2, ...n) . (2) 

1 For the geodesic equations (1) and the formulas applied below, we refer to the standard texts on general relativity. 
In particular, in the following we shall refer to two of the most widespread texts [5, 6]. 
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If U, V are any real functions of the x l (i — 1, 2, ...n), the invariants defined by 

^ U ^ = ^WW S ^ V * (4) 

are called Beltrami, s differential parameters of the first order. Since the equations U = const, 
V = const represent (n — l)-dimensional hypersurfaces in V n , Ai U represents the squared length 
of the gradient of U as well as of a vector orthogonal to the hypersurface U = const ; for the same 
reason, if A(U, V) = 0, then the two hypersurfaces U = const and V = const are orthogonal. 
The Beltrami's theorem states: 
Let us consider the equation 

Ail7 = l, (5) 

the solution of this equation depends on an additive constant and on N — 1 essential constants a, 
[4]. Now, if we know a complete solution of Eq. (5), we can obtain the equations of the geodesies 
from the following theorem [3, p. 299], [4, p. 59]: when a complete solution of Eq. (5) is known, the 
equations of the geodesies are given by 

< 6 > 

where Pi are arbitrary constants, and the geodesic arclength is given by the value of U . 

We demonstrate the theorem following Beltrami's article [2]: 
Let us return to differential parameters (3), (4). We can associate contravariant components of the 
gradient with the covariant ones 

U' h = g ht U,i . (7) 

Since (2) and (3) are reciprocal quadratic forms (since the coefficients in (3) are reciprocal to the 
ones in (2)), because of their invariance we can write 

ds 2 = Ai U dk 2 , (8) 

where dk is an infinitesimal scalar quantity introduced to keep the differential homogeneity. From 
Eq. (8), we have obviously 

ds 

Moreover, owing to the properties of the reciprocal forms and keeping the differential homogeneity 
expressed by Eq. (9), we can write 

g ih dx h = Uidk, (10) 



which can be rewritten as 
Combining (9) with (11), we get 



dx h =g lh U yi dk = U^ h dk. (11) 
dr h U> h 

(12) 



ds ^A~Tu' 

By multiplying (10) by dx 1 and summing, we shall get 

ds 2 = dUdk, (13) 



3 



where dll is the differential of the function U. 

In an analogous way, if the increments 5x % are the components of a vector indicating another 
direction from that of components dx\ we shall get by multiplying by 5x z and summing 

gihSx* dx h = 5U dk . (14) 

It is evident that, if we take the vector indicated by 5x l tangent to a hyper surf ace U = const , we 
have 5U = and then the vectors, of components dx % and b~x % respectively, turn out to be orthogonal. 
Then the 11 displacement" dx of (10), is orthogonal to the hypersurface U = const. 
Finally, if we eliminate dk from Eqs. (8, 13), we find 

where dU represents the increment due to a variation ds which, for what we have seen, is orthogonal 
to the hypersurface U = const. 

Let us now consider a curve represented in parametric form as a function of a parameter t and 

set 

s = ds/dt and x % = dx 1 / dt . 

From (2), one has 

s 2 = g ih x i x h . (16) 

For this quadratic form we can repeat the preceding arguments and, in particular, in place of (10) 
we shall have 

1 s 2 s 

7 = s—— = g ih x h = U :i k;hri = l,...,n, (17) 



2 d x % dx 1 
where now k is a finite constant of proportionality. 

Furthermore, since s 2 and AiU are two scalar quantities in the same metric, we can write 

s 2 = A 1 Uk 2 , (18) 

from which 

k = J_ . (19) 



By eliminating k from (17) and (19) we obtain 

ds Ui 



dx i y/AJJ 



(20) 



The condition for a curve to be a geodesic, as we know, in Lagrangian form (where C = s is the 
Lagrangian) is given by the n equations 

Let us assume that we know n first integrals of these equations 

.n.< J <-) (22) 



4 



and substitute them into s = y gihX 1 x h and ds/dx 1 of (21) which, in this way, will turn out in the 
end to be functions only of the variables x l . If we indicate by (A) the "total" derivatives with 



respect to x % and make explicit the derivatives of the two sides of (21), we get the equations 

/ ds\ d s d s d x r .. 

+ 7T^7n7> (23) 



dx l J dx % dx r dx 
d ( ds\ d d s 



dt \ dx l ) d x r d x % 



x r , (24) 



where the implicit dependence on x l , once one have substituded the first integrals into the expressions 
of s and ds/dx 1 , has been exploited. 

Now, if we differentiate with respect to x 1 the identity s = (ds/dx r ) x r , we get 

. r d ds ds dx r . . 

~ X dx i dx r + dx r dx i ' ^ ' 




From a comparison of (23) and (25) it follows that 

ds d ds 



x r . (26) 



dx i dx i dx r 

These equations and (24) allows us to write the system (21) in the pfaffian form 

( d ds d ds 



x r = 0, (i = l,...,n). (27) 



\dx r dx i dx i dx r / 

The system (27) is satisfied if it is possible to assign a function U (potential) such that 

ds dU 

U r . (28) 



d x r dx r 

From (20) it follows that this condition is satisfied if U is a solution of Eq. 5. Then, if we find a 
solution of the partial differential equation (5), from (12) we obtain the equations of the geodesies 
as functions of s by solving a system of first order ordinary differential equations. If (5) is satisfied, 
we have from (15) that U = s and then U represents the line element along the geodesies. 

On the basis of what we have summarized above, we can finally enunciate the fundamental result 
of Beltrami due to which the second integration can be avoided: If we know a complete solution 
of the partial differential equation (5), we can obtain the geodesic equations by differentiation steps 
alone. 

The demonstration is straightforward. If (5) is satisfied, from (15) it follows that the length of 
the orthogonal segments between two hypersurfaces U = const is the same. According to a theorem 
due to Gauss, the lines which cross orthogonally the fields U = const are geodesic lines and the 
hypersurfaces U = const are said to be geodesically parallel. 

Let us consider a complete solution of (5) which, in addition to an obvious additive constant, will 
contain n — 1 other arbitrary constants cfy. By differentiating (5) with respect to these constants, we 
obtain 

«££ s2 A(tf,g=0,faJ = l n-l, (29) 

which tells us that the hypersurfaces VJ = dU /da\ = const and U = const are orthogonal. If we now 
put 

dU 

V^ = — = /3,, for/ = l,...,n-l, (30) 
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the curves of intersection of these hypersurfaces, orthogonal to U — const, are geodesies and have 
been obtained through differentiation, without being obliged to solve the differential equations (12). 
The demonstration of Beltrami's theorem is complete. 



Obviously the above theorem is particularly useful for applications when we have a complete 
solution of (5) at our disposal or this solution is easily obtained. 

A general case in which the solution of Eq. (5) is obtained in integral form has been obtained by 
Bianchi [3, 4] 2 : 

If the fundamental form (2) can be re-expressed in the generalized Liouville form 



ds 2 = [Xi(an) + X 2 {x 2 ) + ... + X N (x N )} ^ e^dx 1 ) 2 , 

i=i 

where e$ = ±1 and X i is a function of x l alone, a complete integral of Eq. (5) is 

U = C + J2 I \/eJjQ + aij dx\ 

i=l J 

where c and on are constants, the latter being subject to the condition X^=i a i = 0- 
In this case the geodesic equations (30) are immediately given by 



(31) 



(32) 



dll _ 1 r e%dxi 
dai ~2J L. ( Xl + ai ) 



Pi- 



rn 



3 The geodesic equations for the Schwarzschild and Kerr 
metrics 

3.1 The Schwarzschild metric 

We start from the standard form of the so called Schwarzschild metric [8] 



ds 2 '-^df 1 - ■' 2 - ^ ' - 2 °-> ' 



r r — a 

3 



d r 2 - r 2 (d9 2 + sin 2 9 dej) 2 ); a = 2MG 



(34) 



Applied to this metric Eq. (5) gives 

2 



d- 



r — a \ dt 



r — a Idr 
r \dr 



1. 



If we set 



t = A 1 t + A 2 (j) + T 1 (r) + T 2 (9), 
Eq. (35) can be solved by the method of "separation of variables." We obtain 



r — a 



r — a (dr\ 
r \ dr 



dr 2 
~d9 



3 

sin 2 9 



(35) 



(36) 



(37) 



2 In this paper we do not apply this result and obtain the solution of Eq. (5) directly for metrics that are not in the 
Liouville form. 

3 We replace the previous symbol U with r, that recalls the physical meaning of proper time along the geodesies. 
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Since the right hand side is positive, we can set both sides equal to a separation constant A\ to 
obtain 



A- 



r — a ( d T\ 



r — a 



' dr 2 

1e 



dr 
2 



+ 



- 1 



3 

sin 2 9 



A 2 



A 2 



(38) 
(39) 



which then gives 



-I 



^JA\ r 4 - r (r - a)(r 2 + Aj) 
r (r — a) 

JX 2 sin 2 9 - A\ 



T 2 



-I 



sin6> 



dr, 



d6. 



(40) 



(41) 



If in the right hand sides of (40) and (41) we choose the positive sign, the geodesic equations (6) 
become 



3Tr t+A 'I 



r 3 dr 



(r - a)^Alr A - r{r 2 + A 2 3 ){r - a) 



Bi, 



dr 
dAo 



= - A 2 [ - 



d6 



dr 



OA 



-A, 



sm9 J A 2 sin 2 9- A 2 
dr 



= + sin 1 [sinh e cot 9] = B 2 , 



(42) 
(43) 



+ A S 



sin 9 d9 



y/ 'A 2 ^ - r(r 2 + A 2 )(r - a) ~ J \J A 2 sin 2 9 - A\ 

/dr 
, — sin~ 1 [coshe cos^] = B 3 , 

• IA\r A - r(r 2 + Af)(r - a) 



(44) 



where in the explicitly evaluated integrals we have set A% = A coshe, A 2 = A sinhe. 

Eq. (43) is the same as Eq. (4) on p. 645 of Misner et al. [5] which establishes the planar character 
of the orbit. In fact, with a suitable rotation of the angle (0o = B 2 ), we obtain from (43) the result 
sin 4> = — sinh e cot 9 and by substituting this into the defining equations for polar coordinates, one 
finds the orbit in x, y, z space as a function of the parameter 9: 



x = r{9) sin 9 cos <p = r{9) sin 9 y 1 — sinh 2 e cot 2 9 , 
y = r{9) sin 6* sin0 = —r(6) sinhe cos#, 
z = r(9) cos9 . 

It immediately follows that y = — sinhe z, i.e., the orbit lies in this plane. We remark that Misner 
et al. obtain this result by applying the Hamilton- Jacobi method. 
We can now calculate the proper time r from the identity 



. dr . dr 
r-A 1 ——-A 2 - 



dA 1 



8An 



A, 



dr 
'dA, 



r - A 1 B 1 - A 2 B 2 - A 3 B 3 



(45) 



obtained by substituting the partial derivatives by the constants Bi. In the left hand side we evaluate 
t by using Eq. (36) in which T\ and r 2 are given by Eqs. (40) and (41). We also evaluate the other 
three terms of the left hand side by using the first identities of Eqs. (42), (43) and (44), respectively. 
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We group the integrals in d r, d 9 and see that the resulting integral in d 9 vanishes. Finally we obtain 
r as a function of r via an elliptic integral: 



/t dv 
i ■ (46) 

\JA\r 4 - r(r 2 + Af)(r - a) 

At this point we must do two things: compare our results with the existing literature and at 
the same time identify the constants Ai, A 2 , A 3 that we have introduced with the relevant physical 
constants of the classic approach. To do this, we set 9 = n/2, which is the a priori value everyone 
assumes. Then from (39), A\ = A\ and the sum of (43) and (44) we get 



/dr 
+ const , (47) 
^ A 2 r 4 - r{r 2 + A 2 ){r - a) 

which can be re-expressed as 



/dv 
. + const . (48) 
r 2 ^A\-{l + Al/r 2 ){l-a/r) 

In the same way (42) can be rewritten in the form 

/dv 
+ const . (49) 
(1 - a/r) JA 2 -(l + A 2 Jr 2 )(l-a/r) 



If we compare our equations (49) and (48) with (101.4) and (101.5) of Landau-Lifshitz [9], which 
are obtained starting from the Hamilton- Jacobi equation, we find that the equations are the same, 
apart from the different symbols and units used in them. In the limit a/r C 1 (as is the case for 
planetary motion), (48) becomes 

r dr 

6 = A 3 . + const . (50) 

J r 2 \J(A\ - 1) - A 2 /r 2 

If we compare (50) with the relevant Newtonian equation (see Boccaletti-Pucacco [10], Eq. (2.14), 
p. 131) 

c dr . . 

+ const , (51) 



/ 



r 2 ^2 (h - V{r)) - c 2 /r 2 



where c is the angular momentum per unit mass, one immediately identifies A 3 with the angular 
momentum c. From (39), one has that in the general case 9 ^ n/2 the constant A 2 represents the 
component of the angular momentum along the polar axis. 

We can also obtain the physical meaning of the constant A 3 and the second Kepler law from (46) 
and (47). Assuming "direct orbits" we have: 

d(p _ dcp dx_ _ As 
dr dr dr r 2 

As for the constant A±, in our units (c = 1, m — 1), A\ — 1 must be identified with twice the kinetic 
energy through the standard relativistic relations E 2 — m 2 c 4 = c 2 p 2 and p 2 = 2mT, where T is 
the classical kinetic energy. Thus in our units A\ = E 2 ; a further detailed investigation leads to 
A x = -E. 
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Regarding the advance of the perihelion of Mercury, it is convenient to start from (48) rewritten 
in terms of the variable u — l/r 




By differentiating this with respect to 0, we obtain 

d 2 u la 3 



' Aj-(l + Alu 2 )(l-au) . (53) 



# 2 2A\ 2 

which since a = 2G M can be rewritten 

d 2 u GM 



+ U= o^2 + n aU , (54) 



+ ^ = ^(1 + 3^). (55) 



A comparison with the classic Binet equation (see [10], Eq. (2.24)) shows that the second term on the 
right hand side represents the relativistic correction which accounts for the advance of the perihelion. 
The reader can find the details of the calculations in the well known text by Bergmann [11]. 

As a final remark, we point out that, returning to the geodesic equations in the form of (42, 43, 
44), we can further obtain an expression for r as a function of more general than Eq. (47). In fact, 
without the condition 9 = 7r/2, we have from Eq. (43), 9 as a function of and then, from Eq. (44), 
r as a function of 0. This leads to 

cotg = - sin ^- B2) (56) 
sinh e 

and 

t f, / % ^ n cosh e sin (0 — Bo) . „ N 

sm[/(r) + B 3 ] = KV 2) , (57) 

ysinlre + sin 2 (0-S 2 ) 

where 

dr 



f(r) 



\jA\r 4 - r(r 2 + AQ(r - a) 



is the elliptic integral in dr of Eq. (44). Eqs. (56) and (57) give the geodesies in the space as functions 
of the parameter 0. 

Note that above we have obtained the result that the orbits are planar without invoking the 
spherical symmetry of the field and without setting 9 = n/2. 



3.2 The Kerr metric 

The method based on Beltrami's theorem that we have applied so far to study geodesic motion in 
the Schwarzschild spacetime can clearly be applied to the Kerr spacetime as well. We start from the 
Kerr metric [6] 

, o p 2 A , „ S 2 / , , 2aMr , \ 2 . , „ p 2 , „ „ , n9 

ds 2 = *-—-dt 2 (d<p ^o—dt) sin 2 9 - !--dr 2 - p 2 d9 2 , (58) 

S 2 p \ S 2 / A 

where p 2 = r 2 + a 2 cos 2 9, A = r 2 + a 2 - 2 Mr, S 2 = (r 2 + a 2 ) 2 - a 2 A sin 2 0; M and a are constants 
that in the Newton limit represent the mass and the angular momentum per unit mass. 
As is well known (see [6], p. 289) the Kerr metric reduces to the Schwarzschild metric when the 
constant a = 0. 
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The contravariant form of the metric tensor is 



9 l3 = -o 



( S 2 /A 2aMr/A \ 

-A 

0-1 

V 2aMr/A (a 2 sin 2 9 - A) /A sin 2 9 J 



(59) 



where the entries are arranged following the sequence dt, dr, d6, d<f>, respectively. Eq. (5) with the 
contravariant components of the metric tensor given by (59) turns out to be 



S 2 (dr\ 2 AaMr { dr\ { dr\ A — a 2 sin 9 i dr 



A I dt 



+ 



A V dt I \<90 



Asin 2 # \dt 



-A 



d\ 2 / r\ \ 2 



dr 



dt 



(60) 



For (58), as for (34), the coefficients of the line element do not depend on t and 0, therefore we can 
seek a solution of (5) in the same form of Schwarzschild equation 



r = A 1 1 + A 2 + n(r) + t 2 {9) . 
By substituting into Eq. (60) we obtain 



A 
and 
1 



1 r i 2 

A x (r 2 + a 2 ) + A 2 a 



sin 2 9 



A ± a sin 2 9 + A 2 



-A 



' dri 
dr 



' dr 2 
~dJ 



(61) 



= r 2 + a 2 cos 2 9 , (62) 



Ax (r 2 + a 2 ) + A 2 a 



sin 2 # 



Ax a sin 2 9 + A 2 



+ 



dr 2 
~d9 



+ a 2 cos 2 9 . (63) 



The left hand side is a function of r, the right hand side is a positive function of 9; then introducing 
a new separation constant A\ we obtain the solutions 



Tx 



I 



[Ax (r 2 + a 2 ) + A 2 a] 2 - A(r 2 + A 2 ) 
A 



T2 



I 



{Al - a 2 cos 2 9) sin 2 9 -(Ax a sin 2 9 + A 2 ) 2 
sin 9 



dr, 



d9. 



(64) 
(65) 



If we set 

R(r) = [Ax(r 2 + a 2 )+A 2 a} 2 -A(r 2 + Al), 6(0) = (A 2 3 -a 2 cos 2 0) sin 2 9- (Ax a sin 2 9 + A 2 f , (66) 
Eq. (61) becomes 



d0. 



A - J sin* - (67 > 
The equations for the geodesies can now be obtained by the standard procedure of Eq. (6) 
following from Beltrami's theorem. We get 



dr 
dA~x 



t + 



I 



(r 2 + a 2 )[Ax(r 2 + a 2 ) + A 2 a] ^ r (Ax a sin 2 9 + A 2 )asm9 



dr , r A x (r 2 + a 2 ) + A 2 a , r (Ax a sin 2 9 + A 2 

— — = <p + / a — = dr - / 

dA 2 r J Av^R J 



sin^Ve 



d9 = Bx, (68) 
d9 = B 2 , (69) 
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= A, 



+ 



sin 9 



d9 



B 



3 • 



(70) 



Moreover, we can calculate the proper time with the same procedure we have used for the Schwarzschild 
metric. We start from the expression (45) and obtain 



Now the integral in dd does not vanish, while from (70) we can find 9 as a function of r and then find 
r(r). In any case from Eqs. (69, 70, and 71) we can obtain the relations analogous to those obtained 
for the Schwarzschild metric. 

4 Conclusions 

It turns out, as Beltrami himself was the first to point out, that Beltrami's method is formally 
analogous to the method of integration of the Hamilton- Jacobi equation. On the other hand this is 
even more evident if we consider that the differential parameter A 1 r is formally analogous to the 
expression 



which represents the Hamiltonian and then Ai r = 1 is equivalent to the Hamilton- Jacobi equation. 

However, in spite of this formal analogy, by using Beltrami's theorem we remain in a geometric 
context to obtain the geodesies, i.e., the orbits of a test particle in the gravitational field, without 
being obliged to resort to concepts copied from classical mechanics. In fact it must be remarked that 
very often one uses too freely in general relativity procedures which obtain a true legitimacy only for 
r — > oo or weak gravitational fields. The fact that, a posteriori, things turn out to be correct in the 
nonrelativistic limit does not always remove the ambiguity from certain formulations. 

Finally we note that the "geometrical integration" here described allows us to recover the classical 
conservation laws instead of introducing them a priori. 
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